#include "orbit.hh"

/*
          _______________
 R(t) = \/x(t)^2 + y(t)^2
 
             x(t)
 x"(t) = -k ------
            R(t)^3
 
             [y(t)^2-2x(t)^2]x'(t)-3x(t)y(t)y'(t)
 x'"(t) = -k ------------------------------------
                           R(t)^5
 
 k (((3 x(t)^2)/(x(t)^2+y(t)^2)^(5/2)-1/(x(t)^2+y(t)^2)^(3/2)) x''(t)+y(t) (6/(x(t)^2+y(t)^2)^(5/2)-(30 x(t)^2)/(x(t)^2+y(t)^2)^(7/2)) x'(t) y'(t)+x(t) ((9/(x(t)^2+y(t)^2)^(5/2)-(15 x(t)^2)/(x(t)^2+y(t)^2)^(7/2)) x'(t)^2+(3 (x(t)^2 (y(t) y''(t)+y'(t)^2)-y(t)^2 (4 y'(t)^2-y(t) y''(t))))/(x(t)^2+y(t)^2)^(7/2)))
 
              3x(t)^2 - R(t)^2         6R(t)^2 - 30x(t)^2
 x""(t) = k ( ---------------- x"(t) + ------------------ x'(t)y'(t)y(t) + 
                   R(t)^5                    R(t)^7    
 
 */